/*=============================================================================
*
*   Copyright (C) 2020 All rights reserved.
*
*   Filename:		quadrangle.hpp
*
*   Author: Wang Zhecheng - wangzhecheng@yeah.net
*
*   Date:  2020-01-17 16:54
*
*   Last Editors: Wang Zhecheng - wangzhecheng@yeah.net
*
*   Last modified:	2020-08-04 01:24
*
*   Description:
*
=============================================================================*/
#ifndef I_QUADRANGLE_H
#define I_QUADRANGLE_H

#include "polygon_interface.hpp"
#include "quadrangle_shape_function.hpp"

namespace geoxfem {
namespace geometry {

template <typename DimensionType, typename degree>
struct Quadrangle : PolygonInterface<DimensionType, degree> {
  DEFINE_KEY(Quadrangle)

  static constexpr size_t dimension = traits::dimension<DimensionType>::value;

  using Shape      = ShapeInterface<DimensionType, degree>;
  using Polygon    = PolygonInterface<DimensionType, degree>;
  using point_type = typename traits::point_type<DimensionType>::type;
  using Gauss      = model::GaussPoint<DimensionType>;

  /*constructors and destructor{{{*/

  /**
   * @brief constructor.
   *
   * @param _listPoint
   * @param _pMesh
   */
  Quadrangle(const std::vector<point_type>& _listPoint) : Polygon{_listPoint} {}
  Quadrangle()          = default;
  virtual ~Quadrangle() = default;
  /*}}}*/
  /*misc{{{*/

  /**
   * @brief getNumberofPoint
   *
   * @return
   */
  size_t getNumberOfPoint() const final override {
    return degree::value * 4;
  }
  /**
   * @brief get Class Name
   *
   * @return
   */
  std::string getClassName() const final override {
    return std::string("Quadrangle") + Polygon::poly_degree[degree::value]
           + std::to_string(dimension) + "D";
  }
  /*}}}*/
  /*shape function using reference coordinates{{{*/

  /** calculate standard shape function using local coordinates
   *
   * @param lcoord
   * @return
   */
  VectorXR getStdNByLocal(const VectorXR& lcoord) const final override {
    return quadrangle_shape_function<dimension, degree::value>::getStdNByLocal(
      lcoord, *this);
  }
  /** calculate standard shape function derivation dNdlx using local
   * coordinates
   *
   * @param lcoord
   * @return
   */
  MatrixXR getStddNdlxByLocal(const VectorXR& lcoord) const final override {
    return quadrangle_shape_function<dimension,
                                     degree::value>::getStddNdlxByLocal(lcoord,
                                                                        *this);
  }
  /*}}}*/
  /*shape function using global coordinates{{{*/

  /** also you should get shape function N with global coordinate gcoord
   *
   * @param gcoord
   * @return
   */
  VectorXR getStdNByGlobal(const VectorXR& gcoord) const final override {
    return quadrangle_shape_function<dimension, degree::value>::getStdNByGlobal(
      gcoord, *this);
  }
  /** then global coordinate can be yield by function above
   *
   * @param gcoord
   * @return
   */
  VectorXR global2local(const VectorXR& gcoord) const final override {
    return quadrangle_shape_function<dimension, degree::value>::global2local(
      gcoord, *this);
  }
  /*}}}*/
  /*gauss point related{{{*/

  void getStdGaussPoints(std::vector<Gauss>& gp_list) final override {
    size_t   num_gp = model::i_default_gauss_num::set_get().at("quadrangle");
    MatrixXR lcoord(dimension, num_gp);
    VectorXR lweight(num_gp);

    /// form coordinates and weight for gauss points
    switch(num_gp) {
    case 1:
      lcoord << 0, 0;
      lweight << 4;
      break;
    case 4:
      lcoord << 1. / root3, -1. / root3, -1. / root3, 1. / root3,  //
        1. / root3, 1. / root3, -1. / root3, -1. / root3;
      lweight << 1, 1, 1, 1;
      break;
    case 9:
      lcoord << -root3 / root5, 0, root3 / root5, -root3 / root5, 0,
        root3 / root5, -root3 / root5, 0, root3 / root5,  //
        root3 / root5, root3 / root5, root3 / root5, 0, 0, 0, -root3 / root5,
        -root3 / root5, -root3 / root5;
      lweight << 0.31, 0.49, 0.31, 0.49, 0.79, 0.49, 0.31, 0.49, 0.31;
      break;
    default:
      break;
    }
    /* create gauss point list */
    for(size_t i = 0; i != num_gp; ++i) {
      VectorXR gcoord = this->local2global(lcoord.col(i));
      _real    gweight =
        lweight[i] * this->getJacobByLocal(lcoord.col(i)).determinant();
      gp_list.emplace_back(Gauss{gcoord, gweight});
    }
    return;
  }

  /*}}}*/
  /*io related {{{*/
  /**
   * @brief getVTKTypeID
   *
   * @return
   */
  size_t getVTKTypeID() const final override {
    return getVTKTypeIDImpl<degree>();
  }
  template <typename d>
  size_t getVTKTypeIDImpl() const {
    return 53;
  }
  template <>
  size_t getVTKTypeIDImpl<mpl1>() const {
    return 9;
  }
  template <>
  size_t getVTKTypeIDImpl<mpl2>() const {
    return 23;
  }
  std::string getVTKPointOrder() const final override {
    return get_vtk_point_order_impl<degree>();
  }
  template <typename d>
  std::string get_vtk_point_order_impl() const {
    throw "not implemented yet";
  }
  template <>
  std::string get_vtk_point_order_impl<mpl1>() const {
    return Shape::getVTKPointOrder();
  }
  template <>
  std::string get_vtk_point_order_impl<mpl2>() const {
    /* Create an output string stream */
    std::ostringstream streamObj;
    /* Add double to stream */
    streamObj << this->getPoint(0).getIndex() << " "  //
              << this->getPoint(6).getIndex() << " "  //
              << this->getPoint(4).getIndex() << " "  //
              << this->getPoint(2).getIndex() << " "  //
              << this->getPoint(7).getIndex() << " "  //
              << this->getPoint(5).getIndex() << " "  //
              << this->getPoint(3).getIndex() << " "  //
              << this->getPoint(1).getIndex();
    /* Get string from output string stream */
    return streamObj.str();
  }

  /*}}}*/
};

}  // namespace geometry
}  // namespace geoxfem

GEOXFEM_REGISTER_TEMPLATIZED(geoxfem::geometry::Quadrangle);

#endif /* I_QUADRANGLE_H */
